#######################
# This creates Figures and tables for Supplementary section 6, comparing results with different treatment of NA values
#######################

#Initial Results, NA's assumed at midpoint

#OLS regression on ATE for question experiment 1
S6A<-svyglm(tf~consequences,design=nd, family = gaussian(identity))
S6B<-svyglm(tf~consequences+gender+as.factor(racethnicity)+age,design=nd, family = gaussian(identity))

#OLS regression on ATE for question experiment 2
S6G<-svyglm(itf~consequences,design=nd, family = gaussian(identity))
S6H<-svyglm(itf~consequences+gender+as.factor(racethnicity)+age,design=nd, family = gaussian(identity))

#Computing ATE's, treating NA's as NA

#Load Data
nNA<-read_dta("n.dta")

#Re-code Force required to preserve trump for NA as NA, not 0
nNA$tf<-ifelse(nNA$uchi2c %in% c(1,2),1,0)
nNA$tf<-ifelse(nNA$uchi2c %in% c(77,98),NA,nNA$tf)

#Re-code Personally would use force to preserve Trump for NA as NA, not 0
nNA$itf<-ifelse(nNA$uchi2e %in% c(1,2),1,0)
nNA$itf<-ifelse(nNA$uchi2c %in% c(77,98),NA,nNA$itf)

##re-code variable (p_over) which assigns question wording in experiments
#doing this to make it easier for interpretation
nNA$consequences<-ifelse(nNA$p_over==1,1,0)

#reset weighting
ndNA <- svydesign(id      = ~vpsu32,
                         weights = ~weight_uchi,
                         strata = ~vstrat32,
                         nest=TRUE,
                         survey.lonely.psu = "adjust",
                         data    = nNA)

#OLS regression on ATE for question experiment 1, treating no answer as NA
S6C<-svyglm(tf~consequences,design=ndNA, family = gaussian(identity))
S6D<-svyglm(tf~consequences+gender+as.factor(racethnicity)+age,design=ndNA, family = gaussian(identity))

#OLS regression on ATE for question experiment 2, treating no answer as NA
S6J<-svyglm(itf~consequences,design=ndNA, family = gaussian(identity))
S6K<-svyglm(itf~consequences+gender+as.factor(racethnicity)+age,design=ndNA, family = gaussian(identity))

##Hard bounds - treat all NA's as strong support

#Load Data - This assumes you have loaded all replication materials as one
#If not, you will need to change the file path
nNASS<-read_dta("n.dta")

#Code non-response as agree force required to preserve Trump
nNASS$tf<-ifelse(nNASS$uchi2c %in% c(1,2,77,98),1,0)

#Code non-response as agree I would use force to restore Trump
nNASS$itf<-ifelse(nNASS$uchi2e %in% c(1,2,77,98),1,0)

##re-code variable (p_over) which assigns question wording in experiments
#doing this to make it easier for interpretation
nNASS$consequences<-ifelse(nNASS$p_over==1,1,0)

#re-weight results
ndNASS <- svydesign(id      = ~vpsu32,
                           weights = ~weight_uchi,
                           strata = ~vstrat32,
                           nest=TRUE,
                           survey.lonely.psu = "adjust",
                           data    = nNASS)

#OLS regression on ATE for question experiment 1, treating no answer as Support for Question
S6E<-svyglm(tf~consequences,design=ndNASS, family = gaussian(identity))
S6F<-svyglm(tf~consequences+gender+as.factor(racethnicity)+age,design=ndNASS, family = gaussian(identity))

#OLS regression on ATE for question experiment 2, treating no answer as Support for Question
S6L<-svyglm(itf~consequences,design=ndNASS, family = gaussian(identity))
S6M<-svyglm(itf~consequences+gender+as.factor(racethnicity)+age,design=ndNASS, family = gaussian(identity))

#summarize results for all models on support force to restore trump (variable $tf)
SS6A<-summary(S6A)
SS6B<-summary(S6B)
SS6C<-summary(S6C)
SS6D<-summary(S6D)
SS6E<-summary(S6E)
SS6F<-summary(S6F)

#condense summaries for each model on support force to restore trump (variable $tf)
TS6AC1<-round(rbind(SS6A$coefficients[2,1],SS6A$coefficients[2,2],SS6A$coefficients[2,4]),3)
TS6AC2<-round(rbind(SS6B$coefficients[2,1],SS6B$coefficients[2,2],SS6B$coefficients[2,4]),3)
TS6AC3<-round(rbind(SS6C$coefficients[2,1],SS6C$coefficients[2,2],SS6C$coefficients[2,4]),3)
TS6AC4<-round(rbind(SS6D$coefficients[2,1],SS6D$coefficients[2,2],SS6D$coefficients[2,4]),3)
TS6AC5<-round(rbind(SS6E$coefficients[2,1],SS6E$coefficients[2,2],SS6E$coefficients[2,4]),3)
TS6AC6<-round(rbind(SS6F$coefficients[2,1],SS6F$coefficients[2,2],SS6F$coefficients[2,4]),3)

#bind condensed summaries together
TS6AR24<-cbind(TS6AC1,TS6AC2,TS6AC3,TS6AC4,TS6AC5,TS6AC6)

#create a row stating dataset
TS6AR1<-c("NA as midpoint","NA as midpoint","NA as NA","NA as NA","NA as Strong Support","NA as Strong Support")

#create a row stating controls
TS6AR5<-c("None","Gender,Race,Age","None","Gender,Race,Age","None","Gender,Race,Age")

#Bind these together
TS6A<-rbind(TS6AR1,TS6AR24,TS6AR5)

#give descriptive row names
row.names(TS6A)<-c("NA treatment","Estimated ATE","Std Err","p-value","controls")

#give descriptive column names
colnames(TS6A)<-c("OLS Model 1","OLS Model 2","OLS Model 3","OLS Model 4","OLS Model 5", "OLS Model 6")

#write out results
write.csv(TS6A,"../Supplementary Results/TableS6A.csv")

#summarize results for all models on support force to restore trump (variable $tf)
SS6G<-summary(S6G)
SS6H<-summary(S6H)
SS6J<-summary(S6J)
SS6K<-summary(S6K)
SS6L<-summary(S6L)
SS6M<-summary(S6M)

#condense summaries for each model on support force to restore trump (variable $tf)
TS6BC1<-round(rbind(SS6G$coefficients[2,1],SS6G$coefficients[2,2],SS6G$coefficients[2,4]),3)
TS6BC2<-round(rbind(SS6H$coefficients[2,1],SS6H$coefficients[2,2],SS6H$coefficients[2,4]),4)
TS6BC3<-round(rbind(SS6J$coefficients[2,1],SS6J$coefficients[2,2],SS6J$coefficients[2,4]),3)
TS6BC4<-round(rbind(SS6K$coefficients[2,1],SS6K$coefficients[2,2],SS6K$coefficients[2,4]),4)
TS6BC5<-round(rbind(SS6L$coefficients[2,1],SS6L$coefficients[2,2],SS6L$coefficients[2,4]),3)
TS6BC6<-round(rbind(SS6M$coefficients[2,1],SS6M$coefficients[2,2],SS6M$coefficients[2,4]),4)

#bind condensed summaries together
TS6BR24<-cbind(TS6BC1,TS6BC2,TS6BC3,TS6BC4,TS6BC5,TS6BC6)

#create a row stating dataset
TS6BR1<-c("NA as midpoint","NA as midpoint","NA as NA","NA as NA","NA as Strong Support","NA as Strong Support")

#create a row stating controls
TS6BR5<-c("None","Gender,Race,Age","None","Gender,Race,Age","None","Gender,Race,Age")

#Bind these together
TS6B<-rbind(TS6BR1,TS6BR24,TS6BR5)

#give descriptive row names
row.names(TS6B)<-c("NA treatment","Estimated ATE","Std Err","p-value","controls")

#give descriptive column names
colnames(TS6B)<-c("OLS Model 1","OLS Model 2","OLS Model 3","OLS Model 4","OLS Model 5", "OLS Model 6")

#write out results
write.csv(TS6B,"../Supplementary Results/TableS6B.csv")



#write out results - stargazer method
stargazer(S6A,S6B,S6C,S6D,S6E,S6F,omit=c("gender","age","racethnicity","Constant"),omit.stat=c("ll","aic"),type="text",out="../Supplementary Results/TableS6A-Stargazer.html")

stargazer(S6G,S6H,S6J,S6K,S6L,S6M,omit=c("gender","age","racethnicity","Constant"),omit.stat=c("ll","aic"),type="text",out="../Supplementary Results/TableS6B-Stargazer.html")